// to calculate adjustment for incidental parameters bias

version 18
clear all

cd "Replication Package"

do "analysis_programs.do"
set scheme burd

use "vignette_experiment.dta" , clear

qui{
capture program drop mylogit_mle
program define mylogit_mle
args lnf Xb
	qui replace `lnf' = ln(invlogit(`Xb')*(1-$ML_y2) + (1-invlogit(`Xb'))*($ML_y2)) if $ML_y1==1
	qui replace `lnf' =  ln(invlogit(`Xb')*($ML_y2) + (1-invlogit(`Xb'))*(1-$ML_y2)) if $ML_y1==0
end	

// drop if selfseg == 7

global sample inrange(pair,11,19) & !everinattentive 
global outregopt landscape plain coljust(lc) varlabels ///
  stats(b se p) nostar sdec(3) // starloc(1) starlevels(10 5 1)
global outregfinal landscape plain coljust(lc) varlabels ///
  stats(b se p) nostar 
global expvar ib0.tech c.comm 1.commabovecurrent ///
  0.opt 1.displayedsecond   

cmrologit pref $expvar ib0.selfseghigh#1.tech ///
   ib0.selfseghigh#c.comm if $sample , cl(id)

global tech_orig = _b[1.tech]
global comm_orig = _b[c.comm]
global commabovecurrent_orig = _b[1.commabovecurrent]
global opt_unemploy_orig = _b[0.opt]
global displayedsecond_orig = _b[1.displayedsecond]
global selfseghigh_tech_orig = _b[1.selfseghigh#1.tech]
global selfseghigh_comm_orig = _b[1.selfseghigh#c.comm]

local subj : display %6.0f  e(N_clust)
local ll: display %9.1f e(ll)    
global subj_orig = `subj'
global ll_orig = `ll'

domarginswithcurrentcomm if opt!=0 & e(sample) , dydx(1.tech) expression( exp(xb()) / ( 1 + exp(xb()) ) * 100 ) at1(selfseghigh=0) at2(selfseghigh=1) statname(ls) fmt(%9.1f)
matrix A0 = r(table)
global ls0_orig = A0[1,3] 
global ls0p_orig = A0[4,3]
global ls0lc_orig = A0[5,3]
global ls0uc_orig = A0[6,3]
global ls1_orig = A0[1,4] 
global ls1p_orig = A0[4,4]
global ls1lc_orig = A0[5,4]
global ls1uc_orig = A0[6,4]

nlcom cwd0 : -_b[1.tech] / _b[c.comm]
matrix A0 = r(table)
global cwd0_orig = A0[1,1] 
global cwd0p_orig = A0[4,1]
global cwd0lc_orig = A0[5,1]
global cwd0uc_orig = A0[6,1]

nlcom cwd1 : - ( _b[1.tech] + _b[1.selfseghigh#1.tech] ) / ///
   ( _b[c.comm] + _b[1.selfseghigh#c.comm] )
matrix A1 = r(table)
global cwd1_orig = A1[1,1] 
global cwd1p_orig = A1[4,1]  
global cwd1lc_orig = A1[5,1]
global cwd1uc_orig = A1[6,1]

*********************************************
**#bootstrap **
*********************************************
local bsnum = 5000
clear
set obs `bsnum'

gen b_tech = .
gen b_comm = .
gen b_commabovecurrent = .
gen b_opt_unemploy = .
gen b_displayedsecond = .
gen b_selfseghigh_tech = .
gen b_selfseghigh_comm = .
gen b_subj = .
gen b_ll = .
gen b_ls0 = .
gen b_ls0p = .
gen b_ls0lc = .
gen b_ls0uc = .
gen b_ls1 = .
gen b_ls1p = .
gen b_ls1lc = .
gen b_ls1uc = .
gen b_cwd0 = .
gen b_cwd0p = .
gen b_cwd0lc = .
gen b_cwd0uc = .
gen b_cwd1 = .
gen b_cwd1p = .
gen b_cwd1lc = .
gen b_cwd1uc = .

local i = 1
local j = 1
while `i' <= `bsnum' {
    local seed = 20250630 + `j'
    set seed `seed'
    preserve
    use "vignette_experiment.dta" , clear
    bsample, cluster(id)
    capture {
        quietly: cmrologit pref $expvar ib0.selfseghigh#1.tech ///
           ib0.selfseghigh#c.comm if $sample , cl(id)
        local tech_coef = _b[1.tech]
        local comm_coef = _b[c.comm]
        local commabovecurrent_coef = _b[1.commabovecurrent]
        local opt_unemploy_coef = _b[0.opt]
        local displayedsecond_coef = _b[1.displayedsecond]
        local selfseghigh_tech_coef = _b[1.selfseghigh#1.tech]
        local selfseghigh_comm_coef = _b[1.selfseghigh#c.comm]
        local subj : display %6.0f  e(N_clust)
        local ll: display %9.1f e(ll)
        
        domarginswithcurrentcomm if opt!=0 & e(sample) , dydx(1.tech) expression( exp(xb()) / ( 1 + exp(xb()) ) * 100 ) at1(selfseghigh=0) at2(selfseghigh=1) statname(ls) fmt(%9.1f)
        matrix A0 = r(table)
        local ls0 = A0[1,3] 
        local ls0p = A0[4,3]
        local ls0lc = A0[5,3]
        local ls0uc = A0[6,3]
        local ls1 = A0[1,4] 
        local ls1p = A0[4,4]
        local ls1lc = A0[5,4]
        local ls1uc = A0[6,4]
        
        nlcom cwd0 : -_b[1.tech] / _b[c.comm]
        matrix A0 = r(table)
        local cwd0 = A0[1,1] 
        local cwd0p = A0[4,1]
        local cwd0lc = A0[5,1]
        local cwd0uc = A0[6,1]
        
        nlcom cwd1 : - ( _b[1.tech] + _b[1.selfseghigh#1.tech] ) / ///
           ( _b[c.comm] + _b[1.selfseghigh#c.comm] )
        matrix A1 = r(table)
        local cwd1 = A1[1,1] 
        local cwd1p = A1[4,1]  
        local cwd1lc = A1[5,1]
        local cwd1uc = A1[6,1]
    }
    restore
    
    if _rc {
        noisily display "Error on iteration `j', trying next seed."
        local j = `j' + 1
        continue
    }
    else {
        noisily display "Success on iteration `j'"
		replace b_tech = `tech_coef' in `i'
		replace b_comm = `comm_coef' in `i'
		replace b_commabovecurrent = `commabovecurrent_coef' in `i'
		replace b_opt_unemploy = `opt_unemploy_coef' in `i'
		replace b_displayedsecond = `displayedsecond_coef' in `i'
		replace b_selfseghigh_tech = `selfseghigh_tech_coef' in `i'
		replace b_selfseghigh_comm = `selfseghigh_comm_coef' in `i'
		replace b_subj = `subj' in `i'
		replace b_ll = `ll' in `i'
		replace b_ls0 = `ls0' in `i'
		replace b_ls0p = `ls0p' in `i'
		replace b_ls0lc = `ls0lc' in `i'
		replace b_ls0uc = `ls0uc' in `i'
		replace b_ls1 = `ls1' in `i'
		replace b_ls1p = `ls1p' in `i'
		replace b_ls1lc = `ls1lc' in `i'
		replace b_ls1uc = `ls1uc' in `i'
		replace b_cwd0 = `cwd0' in `i'
		replace b_cwd0p = `cwd0p' in `i'
		replace b_cwd0lc = `cwd0lc' in `i'
		replace b_cwd0uc = `cwd0uc' in `i'
		replace b_cwd1 = `cwd1' in `i'
		replace b_cwd1p = `cwd1p' in `i'
		replace b_cwd1lc = `cwd1lc' in `i'
		replace b_cwd1uc = `cwd1uc' in `i'
		
		local i = `i' + 1
		local j = `j' + 1
    }
}

compress
save "incdpar_bootstrap.dta", replace

*********************************************
**#bootstrap **
*********************************************

use "incdpar_bootstrap.dta", clear

gen tech_diff = b_tech - $tech_orig
su tech_diff, d
scalar tech_bc = $tech_orig - r(p50)
scalar tech_se = r(sd)

gen comm_diff = b_comm - $comm_orig
su comm_diff, d
scalar comm_bc = $comm_orig - r(p50)
scalar comm_se = r(sd)

gen commabovecurrent_diff = b_commabovecurrent - $commabovecurrent_orig
su commabovecurrent_diff, d
scalar commabovecurrent_bc = $commabovecurrent_orig - r(p50)
scalar commabovecurrent_se = r(sd)

gen opt_unemploy_diff = b_opt_unemploy - $opt_unemploy_orig
su opt_unemploy_diff, d
scalar opt_unemploy_bc = $opt_unemploy_orig - r(p50)
scalar opt_unemploy_se = r(sd)

gen displayedsecond_diff = b_displayedsecond - $displayedsecond_orig
su displayedsecond_diff, d
scalar displayedsecond_bc = $displayedsecond_orig - r(p50)
scalar displayedsecond_se = r(sd)

gen selfseghigh_tech_diff = b_selfseghigh_tech - $selfseghigh_tech_orig
su selfseghigh_tech_diff, d
scalar selfseghigh_tech_bc = $selfseghigh_tech_orig - r(p50)
scalar selfseghigh_tech_se = r(sd)

gen selfseghigh_comm_diff = b_selfseghigh_comm - $selfseghigh_comm_orig
su selfseghigh_comm_diff, d
scalar selfseghigh_comm_bc = $selfseghigh_comm_orig - r(p50)
scalar selfseghigh_comm_se = r(sd)

gen cwd0_diff = b_cwd0 - $cwd0_orig
su cwd0_diff, d
scalar cwd0_bc = $cwd0_orig - r(p50)
scalar cwd0_se = r(sd)

gen cwd1_diff = b_cwd1 - $cwd1_orig
su cwd1_diff, d
scalar cwd1_bc = $cwd1_orig - r(p50)
scalar cwd1_se = r(sd)

gen ls0_diff = b_ls0 - $ls0_orig
su ls0_diff, d
scalar ls0_bc = $ls0_orig - r(p50)
scalar ls0_se = r(sd)

gen ls1_diff = b_ls1 - $ls1_orig
su ls1_diff, d
scalar ls1_bc = $ls1_orig - r(p50)
scalar ls1_se = r(sd)

*********************************************
**#bootstrap **
*********************************************

clear
set obs 11
gen Parameter = ""
gen mlval = .
gen bcval = .
gen se = .
gen ci_l = .
gen ci_u = .

replace Parameter = "1.tech" in 1
replace mlval = $tech_orig in 1
replace bcval = tech_bc in 1
replace se = tech_se in 1
replace ci_l = tech_bc - 1.96*tech_se in 1
replace ci_u = tech_bc + 1.96*tech_se in 1

replace Parameter = "c.comm" in 2
replace mlval = $comm_orig in 2
replace bcval = comm_bc in 2
replace se = comm_se in 2
replace ci_l = comm_bc - 1.96*comm_se in 2
replace ci_u = comm_bc + 1.96*comm_se in 2

replace Parameter = "1.commabovecurrent" in 3
replace mlval = $commabovecurrent_orig in 3
replace bcval = commabovecurrent_bc in 3
replace se = commabovecurrent_se in 3
replace ci_l = commabovecurrent_bc - 1.96*commabovecurrent_se in 3
replace ci_u = commabovecurrent_bc + 1.96*commabovecurrent_se in 3

replace Parameter = "0.opt" in 4
replace mlval = $opt_unemploy_orig in 4
replace bcval = opt_unemploy_bc in 4
replace se = opt_unemploy_se in 4
replace ci_l = opt_unemploy_bc - 1.96*opt_unemploy_se in 4
replace ci_u = opt_unemploy_bc + 1.96*opt_unemploy_se in 4

replace Parameter = "1.displayedsecond" in 5
replace mlval = $displayedsecond_orig in 5
replace bcval = displayedsecond_bc in 5
replace se = displayedsecond_se in 5
replace ci_l = displayedsecond_bc - 1.96*displayedsecond_se in 5
replace ci_u = displayedsecond_bc + 1.96*displayedsecond_se in 5

replace Parameter = "1.selfseghigh#1.tech" in 6
replace mlval = $selfseghigh_tech_orig in 6
replace bcval = selfseghigh_tech_bc in 6
replace se = selfseghigh_tech_se in 6
replace ci_l = selfseghigh_tech_bc - 1.96*selfseghigh_tech_se in 6
replace ci_u = selfseghigh_tech_bc + 1.96*selfseghigh_tech_se in 6

replace Parameter = "1.selfseghigh#c.comm" in 7
replace mlval = $selfseghigh_comm_orig in 7
replace bcval = selfseghigh_comm_bc in 7
replace se = selfseghigh_comm_se in 7
replace ci_l = selfseghigh_comm_bc - 1.96*selfseghigh_comm_se in 7
replace ci_u = selfseghigh_comm_bc + 1.96*selfseghigh_comm_se in 7

replace Parameter = "CWD low skill" in 8
replace mlval = $cwd0_orig in 8
replace bcval = cwd0_bc in 8
replace se = cwd0_se in 8
replace ci_l = cwd0_bc - 1.96*cwd0_se in 8
replace ci_u = cwd0_bc + 1.96*cwd0_se in 8

replace Parameter = "CWD high skill" in 9
replace mlval = $cwd1_orig in 9
replace bcval = cwd1_bc in 9
replace se = cwd1_se in 9
replace ci_l = cwd1_bc - 1.96*cwd1_se in 9
replace ci_u = cwd1_bc + 1.96*cwd1_se in 9

replace Parameter = "Work low skill" in 10
replace mlval = $ls0_orig in 10
replace bcval = ls0_bc in 10
replace se = ls0_se in 10
replace ci_l = ls0_bc - 1.96*ls0_se in 10
replace ci_u = ls0_bc + 1.96*ls0_se in 10

replace Parameter = "Work high skill" in 11
replace mlval = $ls1_orig in 11
replace bcval = ls1_bc in 11
replace se = ls1_se in 11
replace ci_l = ls1_bc - 1.96*ls1_se in 11
replace ci_u = ls1_bc + 1.96*ls1_se in 11
}

gen pval = 2 * (1 - normal(abs(bcval / se)))

format mlval bcval se ci_l ci_u pval %9.3f
compress
save "incdpar_bootstrap_result.dta", replace
